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Abstract: The origin of the e + e _ 511 keV line observed by INTEGRAL remains unclear. 
The rate and morphology of the signal have prompted questions as to whether dark matter 
could play a role. We explore the case of dark matter upscattering in the framework of 
eXciting Dark Matter (XDM), where WIMPs x> interacting through a new dark force, 
scatter into excited states x*> which subsequently emit e + e~ pairs when they de-excite. 
We numerically compute the cross sections for two Yukawa-coupled DM particles upscat- 
tering into excited states, specifically considering variations motivated by recent N-body 
simulations with additional baryonic physics. We find that that I > components of the 
partial-wave decomposition are often significant contributions to the total cross section and 
that for reasonable ranges of parameters dark matter can produce the ~ 10 43 e + /s observed 
by INTEGRAL. 
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1 Introduction 

Over the past several years there has been increasing evidence for a variety of astrophysical 
anomalies. These anomalies have generally taken the form of the presence of a new signal 
of radiation or cosmic rays, beyond what was conventionally expected. They come in 
the form of high energy e + e~ sources, as seen by PAMELA [1] and Fermi [2], microwave 
emission from the galactic center [3, 4], and diffuse gamma rays from a broad (20 — 40°) 
range around the galactic center [3]. All of these anomalies can be related to the presence 
of a new, primary source of high energy (~ 100 GeV) e + e~, which may be attributable to 
a weak-scale dark matter origin, either through annihilation or decay. 

An outlier in the list of astrophysical anomalies is the INTEGRAL 511 keV signal [5, 6] 
(see [7] for a recent discussion). Both bulge and disk-correlated sources are observed. The 
morphology of the bulge component is best modeled by a combination of gaussians with 
widths 3° and 11°. The bulge and disks fluxes are comparable at roug hly 10~ 3 ph cm J s . 
This corresponds to a positron production rate of 10 43 e + /s. While this signal has persisted 
over decades, originally seen in the early 1970's [8-10], the origin of the enormous source 
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of positrons needed to explain it remains elusive. In particular, the 10 43 e + /s concentrated 
in the galactic center region, as well as the highly spherical morphology are a challenge to 
achieve from most galactic sources, which tend to trace the disk. Alternative explanations, 
such as low mass X-ray binaries (LMXBs) [11] could possibly provide a candidate [12], 
although no point source 511 keV emission has yet been observed [13]. Likewise, it has 
been suggested that the transport of the positrons produced in the galactic disk into the 
galactic center could provide the rate [14], although the precise dynamics that achieves this 
and yields such a spherical morphology is unclear. 

At the same time, the possible connection of this signal to dark matter is even less 
obvious, principally because of two facts: first, that the shape of the 511 keV line is 
sufficiently narrow as to constrain the injection energy to be below ~ lOMeV [15]. Second, 
for a weak scale dark matter candidate, the rate is orders of magnitude above what is 
expected from a thermal WIMP annihilation signal. Such ideas have induced people to 
focus on MeV scale dark matter particles [16-18] as an alternative, but a connection to 
more massive theories of dark matter remains appealing. 

A candidate explanation of this was the "eXciting Dark Datter" (XDM) proposal [19]. 
In this proposal, an excited state x* °f the dark matter x 1S postulated with a splitting 
5 ~ 2m e 1 . Dark matter - dark matter collisions mediated by a new dark force produce 
an excitation x ~ > X* followed by the decay x* ~~ ^ X e+e ~ ■ The appealing aspect of this 
idea is that one converts the kinetic energy of a WIMP into positrons. Since this is a 
scattering, rather than annihilation process, the cross section can be much larger, giving 
the possibility of yielding the enormous O(10 43 e + /s) observed in the galactic center. 

Such a proposal is not without problems however. In order to produce the large rates, a 
large cross section a ~ l/c/ 2 was needed, forcing the inclusion of a new GeV-scale mediator, 
(ft, and even then remains challenging [20-22]. Intriguingly, because XDM freezes out by 
annihilating into (ft, the annihilation XX ~^ W1U naturally produce a hard positron signal 
[19, 23]. This simple idea led to the proposal of a "unified" model [24], which simultaneously 
addresses PAMELA/Fermi, WMAP, INTEGRAL and DAMA (through the inelastic dark 
matter scenario [25]). It was argued that in such models (ft can mediate a Sommerfeld 
enhancement of the annihilation [24, 26] 2 to give the rates observed at PAMELA. 

Our focus here is to reconsider more carefully the low-energy positron signal. Because 
the signal involves a non-perturbative process, mediated by multiple (ft exchanges, calculat- 
ing the expected rates can have important subtleties. Moreover, the possible rates depend 
sensitively on both astrophysical and particle physics parameters. By approaching this in 
detail, we hope to understand what the natural expectation for a rate in the galactic center 
would be and what ranges of astrophysical and model parameters can explain the observed 
INTEGRAL signal. 

1.1 Models of XDM and Signals at INTEGRAL 

Models of XDM are simple to construct [19]. The first proposed model involved a dark 
matter as a complex scalar x coupled to a new gauge boson (ft^ of a dark gauge group 

1 For related work, see [20]. 

2 The Sommerfeld enhancement was first discussed in the context of dark matter by [27, 28]. 
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Figure 1. Process responsible for the INTEGRAL signal in the XDM model with a vector mediator. 
(left) Scattering excitation process for XDM is enhanced by multiple (j> exchange (not shown), (right) 
The excited states decays back to the ground state through an offshell 0, producing e + e~ pairs. 

U(l)d, with a small ^ GeV mass. We assume that after U{l)d breaking a small splitting 5 
arises between the real scalar states. 

£ = (D M Xt)*AiXi + \ F U FdilV + ^F AilV + + M 2 X *Xi + M&xiXi + h.c.(l.l) 

One can also replace the scalar easily with a pseudo-Dirac fermion, 

C = i Xi Pxi + \f%F* w + eF^F*" + m 2 ^r + M^iXi + $iXiXi + h.c. (1.2) 

The relic abundance is established by freezing out through annihilations into 0, which 
yields the usual "WIMP miracle" result that for (av) «3x 10~ 26 cm 2 , tth 2 0.1. 

The light force carrier couples to SM states arises through kinetic mixing, as first 
described by [29]. Although this mixing can be quite small when only addressing the 
INTEGRAL signal, for sizeable e, terrestrial experiments (such as fixed target, beam dump 
and searches at low-energy accelerators) can provide limits on these light bosons [30-32]. 
Indeed, there are already important new limits [33-36] and additional proposals for new 
searches [37-41]. 

Even with this setup and a light mediator, it is not obvious that such a model can 
actually achieve such a large rate. We can estimate this by using the Einasto profile [42] 
with parameters set by the A-l run of the Aquarius simulation [43]. In the presence of a 
light mediator, a natural scale for the scattering rate is set by the geometric cross section 
a ~ ft/q 2 - At threshold, q 2 = m x 5, which we take as a lower bound on the natural scale 
of the scattering. Assuming a relative velocity of 2 x 10 _3 c, m x ~ 1 TeV and S ~ 1 MeV, 
one estimates a rate in the inner 2 kpc of 7 x 10 42 e + /s. 

This can be increased, for instance for particles well above threshold, where the mo- 
mentum transfer is low. Our order of magnitude calculation, however, assumes that all 
particles are kinematically capable of scattering, which is an overestimate. For the U(l) 
vector interaction, both WIMPs are excited, requiring 25 of available energy. For this to 
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occur, the characteristic velocity of a 1 TeV WIMP must be ^ 425 km/s. If the velocity 
dispersion is low ~ \/3/2 x 220 km/s - comparable to most estimates of the local value - 
only a small fraction of particles will be kinematically capable of scattering. Indeed, with a 
low and constant velocity dispersion, it is essentially impossible to achieve these high rates 
[20]. Such observations prompted the development of scalar mediated models [19] and 
non-Abelian models [22, 24, 44, 45], as well as models with metastable states [22, 46-48], 
where the threshold is lower. 

However, it would be surprising if the velocity dispersion would stay constant, and a 
number of recent simulations with baryons [49-53] see an increase of the dispersion roughly 
as a power law of the radius as one moves towards the galactic center. If this is true, then 
in the inner 1 kpc the majority of particles would be capable of scattering and the high 
velocities can allow for larger cross sections and scattering rates. Even so, it is not clear 
that such high rates can be achieved, with [22] finding no points in parameter space that 
can achieve these high rates. 

In this paper, we will re-examine this question. In section 2 we will discuss our approach 
to calculating the scattering process, which is consistent with previous approaches using 
a partial wave analysis. In section 3 we convert these partial wave amplitudes into the 
expected rates for INTEGRAL and explore the contributions from each partial wave mode. 
Using Einasto parameters from the Aquarius A-l DM-only simulation [54] we find rates of 
10 -10 42 e + /s. In this section we also explore variations of individual profile parameters 
and find that they can change the rates by a factor of 5-10. In section 4 we consider more 
recent simulations including baryons [52, 53] and find rates of 10 -10 43 e + /s, enough to 
explain the excess. Here we also compare our work with previous work on this matter. 
Finally, in section 5, we discuss connections to other signals and conclude. 

2 Calculational Approach 

Let us consider a two-state system where the states are separated by a mass splitting 5. 
We are interested in 2 — > 2 scattering where two particles enter in the ground state and 
both are upscattered into the excited state. In this scenario, the total splitting between the 
incoming two-particle wavefunction and the outgoing two-particle wavefunction is defined 
to be A. To avoid confusion we will only refer to this total splitting, A. Note that this 
total splitting between the two-particle wavefunctions due to a double excitation, where 
A = 25, is equivalent to a single excitation with A = 5. 

Because the particles are moving non-relativistically, the system is simply governed by 
the Schrodinger equation, which we will solve in the basis of partial waves. We assume the 
particles are attracted by a Yukawa-type force mediated by a particle with mass jtu. We 
will consider XDM-type scenarios where the coupling is off-diagonal. Note that we use 
for the fine structure constant so as not to confuse it with the Einasto profile parameter 
a. For each partial wave mode I, the reduced Schrodinger equation has the form 

v .(»fA + (>JW- E )(*fA (2 ,) 

m x \Xi(x) J \X2(x) J V m x r J \X2{x) J 
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where E is the energy of the two-state system and the potential V is given by 



V= "- ad —\ (2.2) 



Following [24] we restate the Schrodinger equation into the dimensionless parameters 
e v = v/ctd, e<5 = \/ A/m x /ad and = m^/fa^). Using this reparameterization (and 
rescaling r by r — > adm x r) we are left with 

(2.3) 

For boundary conditions, we chose the wave functions to be regular at the origin. 
Remember that for the reduced Schrodinger equation: x( r ) = rR(r), where R(r) is the 
radial part of the solution to the full (spherically symmetric) Schrodinger equation. This 
gives the two conditions xi(0) = and X2(0) = 0. We then impose that X2 is composed of 
purely outgoing spherical waves at infinity. This leaves us with one condition left and to 
set it we simply normalize xi = 1 a t infinity. We are free to do this because we will always 
be concerned with ratios of the wavefunctions. 

We would like to note two things about these expressions. The first is that we can see 
this parameterization is equivalent to that of [21] by noting that T = 1/e?, T = (e^/e|) 2 
and 7] = (e^/e|). The second is that the equations depend only on v, ad and the ratios 
A/m x and m ( k/m x . In this work we will only consider = 1/100 and total splittings of 
1 MeV and 2 MeV. Also, since we are concerned only with thermalized cross sections, the 
velocity will always be integrated over. This leaves us with only two physical parameters 
to scan over: the WIMP mass and the force carrier mass. 

We would also like to reiterate the statement of [55] that much of the parameter space 
is numerically unstable and thus we found it difficult to accurately compute partial wave 
modes higher than I = 7. We found the most stable computational method to be a variant 
of the shooting method called the chasing method [56, 57] with roughly 50 digits of precision 
during the internal computation. We imposed strict error tests on the auxiliary system 
for the unknown boundary value and rejected any data points that failed those tests. In 
section 2.1 we will present the results of these computations as plots of the partial wave 
modes. In section 3.2 we discuss the convergence of our sums over partial waves, but we 
note that if anything this technique underestimates the rates by truncating the sum at 
I = 7. 

2.1 Partial Waves 

We are ultimately interested in the rate of e + e _ pair production, which depends on the 
thermalized scattering cross section. But to compute the cross sections for upscattering, 
we must first compute the partial wave amplitudes, f\. For a given a^, A/m x and m^/m x , 
each fi is a function of e v . We find that values of e v greater than about 0.35 are highly sup- 
pressed, independent of the WIMP characteristics. This constraint comes from assuming a 
maximum escape velocity of 1000 km/s in the center of the galaxy. As discussed in section 
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2.3, the low end of the velocity integral is set by the threshold velocity vth 
Thus we need only compute the fi functions in this window. 
We define the partial wave amplitudes as 



= 2^/A/m. 



^1X2. 

k\xi 



2 



fl 



out 



(2.4) 



2 



where k' is the momentum of the excited state and k is the momentum of the ground 
state. This is the same definition used in [21]. It can be obtained by using the conservation 
of probability flux to define the scattering amplitude. We numerically separate the wave 
functions into incoming and outgoing components by taking Fourier transforms. 

To get a sense of these partial wave amplitudes we present a selection of them as 
functions of e v for various masses and splittings. They have similar shapes to the partial 
wave amplitude functions found by [22]. Any gaps along the functions where the plot 
marker is missing correspond to a point that yielded some sort of numerical error. While 
the partial wave amplitudes we show here have similar shapes to those of [22] , the parameter 
space considered here does not completely overlap with the parameter ranges considered 
there. We will return to this point and its implications in section 4.1. 

2.2 Velocity Profile 

Now with the partial wave amplitudes in hand, we are nearly ready to construct the 
cross sections. But because the partial wave amplitudes are velocity dependent, our cross 
sections will also depend on velocity. We then must thermally average our cross sections 
over the WIMP velocity distribution in the galaxy, in order to yield accurate results. We 
will first then develop the appropriate velocity distribution before we can calculate the 
thermally-averaged cross sections. 

In the rest frame of the galaxy, we assume the WIMPs to have at every point a Maxwell- 
Boltzmann speed distribution peaked around the RMS speed vq. In principle vq could be a 
function of galactic radius. As stated above, the constant vq case doesn't yield the desired 
rates, so we will consider the case where vo(r) = (220 km/s)( 8 [ pc )~^ 4 . In section 4 we 
will consider velocity dispersions specific to a set of DM simulations including baryons as 
well. The WIMP speed distribution is truncated by the escape velocity, v esc . The escape 
velocity is assumed to be a function of galactic radius. 

To estimate the escape velocity profile, we assume the rotational velocity profile is 
fairly flat as a function of radius [58] . The condition for uniform circular motion then gives 



We can break the integral up into an integral from r to Rq and an integral from R@ to oo. 
The second integral is just our local escape velocity. Plugging in for M(r) from equation 



v\ _ GM{r) 




which we can then plug into the escape velocity condition 




(2.6) 



-6- 



rax = 239 GeV 



mx = 368 GeV 




rax = 569 GeV 



rax = 8? 9 GeV 




0.2 0.3 0.4 0.5 0.6 0.7 0. 



rax = 1000 GeV 



rax = 1357 GeV 




Figure 2. Partial wave amplitudes as functions of e v for various m x . m$ = 1 GeV and A = 1 MeV 



2.5 and assuming a local escape velocity of 600km/s [59] we are left with the escape velocity 
distribution 

(2.7) 



2 -27; c 2 ln( ^) +(600km/s) 2 



This allows us to know the escape velocity in terms of the circular velocity v c with a 
reasonable assumption of the total (dark+baryonic) mass distribution function for the 
galaxy. 

We will work in the center of momentum frame and will thus need to transform the two 
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Figure 3. Partial wave amplitudes as functions of e v for various m^. m x = 1 TeV and A = 1 MeV 



particles' three-dimensional velocity distributions into a one-dimensional relative velocity 
distribution. To do this, we first change variables from the two particles' velocities {v{ 
and V2) to the total and relative velocities: v"t = v{ + V2 and v r = v{ — vi- We can then 
integrate out the angular parts and radial part of the total velocity with the conditions 
that v{, V2 < v esc . This leaves us (after proper normalization) with a velocity distribution 
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Figure 4. (left) RMS velocity profile, (right) escape velocity profile, (dot-dashed) v c = 220km/s, 
(solid) v c — 250 km/s. 
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Figure 5. Escape velocity profiles for different local escape velocities, dotted: vi oc — 400 km/s, 
dot-dashed: vi oc = 500 km/s, dashed: vi oc — 600 km/s and solid: vi oc = 700 km/s. 



that is solely a function of v r : 



-4 

e 2v o v. 



e 2^rJ _ e "o 



V( Vr 



v V^v e~^erf(-^-) - 2v, 



which has been normalized by two factors of N, where 

N = 47rv 2 e- v2 /^ dv. 

Jo 



'(2v esc - V r ) 



(2.8) 
(2.9) 



We include for reference in figure 4 plots of the RMS velocity profile and the escape 
velocity profile for two different local circular velocities. In this work we will generally 
assume a local circular velocity, v c , of 250 km/s [58]. We also include in figure 5 the escape 
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velocity profiles for v c = 250 km/s and various local escape velocities. We will consider 
the effects these escape velocity profiles have on the rates in section 3.3.3. We plot the 
profiles up to a radius of r ~ 1.2 kpc, which roughly corresponds to the angular size of 
the INTEGRAL signal. The dashed vertical line in the plots marks the lower limit on our 
integral which we have chosen to be r = 0.075 kpc (to avoid the cusp as r — > 0). 

2.3 Cross Sections 

Armed with the partial wave amplitudes and the relative velocity distribution we are finally 
ready to construct the thermalized cross sections. Standard partial wave scattering theory 
tells us the cross section for a given I mode is given by 



(7,(6,,) = J dfi|(2J + l)Al5(coBI 

n{2l + l)k'\x2out\ 2 



k' 2 k\Xlir 

tt(2Z + 1) 



(2.10) 



mtv, 



X r 

where e v is v/a^ = (v r /2)/a,i and we have used the ground state mass in the last line 
rather than the excited state. The thermalized cross section is then given by the sum of 
these I partial cross sections integrated over the velocity distribution: 

oo „2v 

( av )=J2 "° <ri(^-)v r V(v r )dv r . (2.11) 

1=0 v th " 

Here the lower limit is given by the threshold velocity for inelastic scattering. The threshold 
velocity, Vth = 2../ — , is the minimum relative velocity needed to scatter when there is a 



mass difference and it satisfies 2 x ^m x (^) 2 = A. 
2.4 Einasto Density Profile 

The rate of scattering per volume of two identical WIMPs is given by 

§ = l<(av) (2.12) 

where n x is the WIMP number density and the factor of 1/2 avoids double-counting. In 
general, n x is a function of galactic radius. Recent halo simulations including the effects 
of baryons [52, 53, 60, 61] favor Einasto density profiles so those are the ones we shall 
consider. 

The generic Einasto density profile is given by 

>(r)\ -2 



log 



r 
r-2 



1 



P-2 J a 

3 



(2.13) 



We can eliminate p_2 by assuming a local WIMP density of 0.4 GeV/cm . This leaves us 
with the density profile 

'2ffR Q \ a (r 



p(r) = - exp 
5 



(2.14) 
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Figure 6. Number densities from the Einasto profile using a — 0.17 and r_2 = f5.79. 

Here a (not to be confused with the dark fine structure constant) determines the 
cuspiness of the profile and r_2 is the radius at which the logarithmic slope takes the 
isothermal value. As a baseline, we will assume a = 0.17 and r_2 = 15.79, which correspond 
to the A-l run of the Aquarius simulation [54], though we will we consider variations of 
these parameters in section 3.3. As a reference, we include plots of the WIMP number 
density for both the 100 GeV and 1 TeV WIMPs for these parameters. 

We can then integrate dT/dV over a volume in the galaxy to get the total rate of 
scatterings in that volume. We chose to integrate in a small region in the center of the 
galaxy with radius r c . Remember that both the density profile and the thermalized cross 
section (via the RMS and escape velocities) depend on galactic radius. The final expression 
for the total rate of scattering, T, is then 



Note that due to the cuspiness and uncertainty in the center of the galaxy, we do not 
actually integrate from — we start the integral at r = 0.075 kpc. 



We numerically solve the Schrodinger equation for many points in the m x -m ( f ) parameter 
space to construct the partial wave amplitude functions. We first hold fixed at 1 GeV 
and vary m x from 100 GeV to 5 TeV. Next we hold m x fixed at 1 TeV and vary from 
500 MeV to 5 GeV. We do this for both A = 1 MeV and A = 2 MeV. We then numerically 
integrate the differential scattering rate from the center of the galaxy to 1.2 kpc with the 
Einasto profile parameters a = 0.17 and r_2 = 15.79, assuming a local DM density of 
0.4 GeV/cm 3 . This distance roughly corresponds to the angular width of the INTEGRAL 
signal. For a detailed discussion of what constraints can be inferred by requiring the dark 
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Figure 7. Rates of e + e production using Aq-A-1 profile parameters, disks: A = f MeV, squares: 
A = 2 MeV. 

matter to fit the angular profile of the signal, we would refer the reader to the more detailed 
discussions in [48, 62]. 

We can see from figure 7 that with these parameters, the rates are of the order O(10 41 ) — 
O(10 42 ). We also see there are resonance regions along variations in m x . There are also 
smaller resonances in m^, but generally lower m^'s give higher rates. Of the two resonance 
regions along = 1 GeV, A = 1 MeV, the one peaked around m x « 200 GeV is probably 
heavily dependent on the velocity profile, as we will show in section 3.2. The peak around 
m x ~ 600 GeV is more robust considering possible changes in the halo model. In section 
3.3 we will show that these rates can easily go up an order of magnitude or more by varying 
the profile parameters. 

3.1 Partial Wave Contributions 

To better understand these rates, let us look at the partial wave composition of each 
rate. In figures 8, 9, 10 and 11 we show the absolute and fractional contribution to the 
total rate from each partial wave. We see that higher values of tend to be dominated 
by very low I modes while the lower masses are more uniformly populated. This makes 
intuitive sense if we think of the Compton wavelength of ms as the characteristic radius for 
the angular momentum. Higher means the two WIMPs must get closer to upscatter 
and thus the system has lower angular momentum. The opposite case is true for the m x 
variations. Here the low m x s are dominated by low I modes and the higher masses are 
evenly populated. This is because the lower the m x , the higher the velocity needed to 
overcome the inelastic threshold. Thus the lower-mass WIMPs are all coming from the tail 
of the velocity distribution. This means that when they do upscatter, they are more likely 
to give up all their kinetic energy and the system is left without any angular momentum. 
A higher-mass WIMP at the same velocity will be capable of higher angular-momentum 
processes, as it can remain in a high I state after the transition due to its relatively higher 
residual kinetic energy. 
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Figure 8. Individual partial wave contributions for = 1 GeV and A = 1 MeV varying m x . 
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Figure 9. Individual partial wave contributions for = 1 GeV and A = 2 MeV varying m x . 
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Figure 10. Individual partial wave contributions for m x — 1 TeV and A = 1 MeV varying m^. 



3.2 Convergence of Sums 

A concern we might have is that we may be significantly underestimating the partial wave 
sums with so few I modes. To examine this, we have looked at how these saturate their 
"total" as we increase the sum from I = 5 to I = 7. We find that when summing up to 
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Figure 11. Individual partial wave contributions for m x = 1 TeV and A = 2 MeV varying m^. 
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Figure 12. Differential and total rates as a function of galactic radius, (solid) Percentage of total 
rate, (dot-dashed) Differential rate (in units of kpc^ 1 ) normalized by the total rate. 



I = 5 over 80% of the points we considered are within 20% of their values from summing 
up to I = 7. With / = 6, over 90% of the parameter points are within 20% and almost all 
of the points have reached at least 70% of the I = 7 total. Moreover, from section 3.1 we 
can see that the rates with the most significant high I contributions are either very high 
m x or very low m^. From this we conclude that in the most relevant regions of parameter 
space, the I > 7 modes do not give significant contributions, and even where they do, we 
expect our results to be correct to 0(1). 

Because of the cuspiness in the profiles, another way the rates might be deceiving 
us is if they reach their total value in the first few hundred parsecs. For example, this 
could imply low mass WIMPs were garnering all their scatterings from questionably high 
velocities in the very center of the galaxy. In figure 12 we plot the fraction of the total rate 
achieved as a function of galactic radius. We see that for the 1 TeV WIMP with a splitting 
of 2 MeV the rate has reasonable contributions from all parts of the integral. The 154 GeV 
WIMP with a 1 MeV splitting on the other hand picks up the majority of its rate in the 
first 300 pc. The accuracy of this rate relies upon the precise details of the RMS velocity 
profile, the escape velocity profile and the Einasto density profile in the very inner region 
of the galaxy. 
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Another way to see this is to ask: at a given radius and for a given A/m x , what 
fraction of the velocity profile kinematically allows upscattering? We can plot this fraction 
as a function of radius to see if only the tail is contributing or if a significant portion of 
the WIMPs are contributing. We see in figure 13 that, at best, only 4 — 6% of the 154 GeV 
WIMPs are upscattering while 30 — 60% of the 1 TeV WIMPs can upscatter. This reaffirms 
our conclusions in section 3.1 that the low mass WIMPs (with only low I contributions) 
are sampling from the tail, while higher mass WIMPs sample a broader range of particles. 
We expect that our results for higher mass WIMPs should be fairly robust. 

mx = 154 GeV A = 1 MeV m = i TeV A = 2 MeV 
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Figure 13. Fraction of particles in the halo which are kinematically accessible for upscattering. 
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3.3 Profile Variations 

For all of the rates shown so far the Aquarius A-l profile values of a = 0.17 and r_2 = 15.79 
have been used. But since there is significant uncertainty in the profile parameters, let us 
consider the effects on our rates from variations of these parameters. In the following 
sections we will look at changes in the two Einasto parameters a and r_2 as well as 
variations in the local escape velocity. 

3.3.1 Variations of a 

In figure 14 we plot the effects on the rates of varying the Einasto profile parameter a. 
We use as an example the scenario with = 1 GeV and A = 1 MeV for various m x . 
We vary a from 0.05 to 0.2 in steps of 0.05 while fixing r_2 = 15.79, v c = 250 km/s and 
vioc = 600 km/s. We see that a variation of a over this range can change the rates by an 
order of magnitude or more with the rate increasing as a decreases. 

3.3.2 Variations of r 2 

In figure 15 we plot the effects on the rates of varying the Einasto profile parameter r_2- 
We use the same example scenario of = 1 GeV and A = 1 MeV while varying m x . We 
fix a = 0.17, v c = 250 km/s and v\ oc = 600 km/s and then vary r_2 from 12 to 21 in steps 
of 3. We see that at best variations in r_2 can change the rates by about a factor of 5, 
with the lower r_2 values giving higher rates. 
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Figure 14. Effects on rates of varying a. dotted: a = 0.05, dot-dashed: a — 0.1, dashed: a = 0.15 
and solid: a = 0.2 
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Figure 15. Effects on rates of varying r_2- dotted: r_2 = 12, dot-dashed: r_2 = 15, dashed: 
?"_2 = 18 and solid: r_2 = 21 

3.3.3 Variations of Local Escape Velocity 

In figure 3.3.3 we plot the effects on the rates of varying the local escape velocity (which 
in turn varies the escape velocity in the center of the galaxy). Again, we use the same 
example scenario of rrij, = 1 GeV and A = 1 MeV while varying m x . Here we fix a = 0.17, 
r_2 = 15.79 and v c = 250 km/s. We plot rates for local escape velocities of 400 km/s, 
500km/s, 600 km/s and 700 km/s. As expected the rates go up for higher escape velocities, 
but the effect is mainly in the lower mass WIMPs. From 400 km/s to 700 km/s we see an 
overall enhancement of about a factor of 5 for mx = 500 GeV and lower masses see an order 
of magnitude or more. In light of figure 12, these low m x enhancements are probably due 
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to contributions in the innermost region of the galaxy. It is interesting to note that rates 
for mx > 1 TeV are roughly independent of the local escape velocity. This is due to the 
fact that the higher mass WIMPs have high percentages of kinematically allowed pairs. 

Variation of local escape velocity 




5000 



m x (GeV) 

Figure 16. Effects on rates of varying the local escape velocity, dotted: v\ oc = 400km/s, dot-dashed: 
vice — 500km/s, dashed: vi oc = 600km/s and solid: vi oc — 700km/s. 



4 Simulations with Baryons 

Dark matter only simulations can probe high resolutions and short distances, but since we 
are most interested in the galactic center, where baryonic physics is important, the most 
reliable thing would be to employ simulations that also include baryonic physics. Recently, 
[53] re-simulated many of the original, DM-only Aquarius runs while including the effects 
of baryons. They find that in the inner, baryon-dominated regions the halos become 
more concentrated. While the formation history plays a significant role in determining the 
characteristics of a galaxy's DM halo, the presence of baryons in the simulations seems 
to significantly increase the DM densities in the inner galactic region and this in turn 
provides significant increases to the e + e~ production rates. While the local DM densities 
and local velocity dispersions are similar to those used throughout this paper, their velocity 
dispersions have a much weaker dependence on radius. In table 1 we give for reference 
various values from these re-simulated scenarios. 

We can then use the profile parameters from these simulations to calculate pair pro- 
duction rates. We present variations over and m x for both A = IMeV and A = 2MeV. 
Now though a, r_2, p~2 and the velocity dispersion are all quantities determined from each 
individual simulation (see table 1). As shown in figure 17 the DM densities in the inner 
part of the galaxy are generally higher than the A-l DM-only simulation while their local 
densities can be higher or lower. For most runs these increases in density give increases of 
roughly 3-40 to the rates, as shown in figures 18 and 19. 
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Table 1. Profile values from the DM plus baryons simulations. 
DM energy densities (inner galactic region) DM f ner gy densities ( loca f) 
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Figure 17. Densities for a 100 GeV WIMP in Aquarius runs re-simulated with baryons (see [53]). 
solid (from top to bottom): Aq-A-5, Aq-E-5, Aq-C-5, Aq-D-5, Aq-F-5, Aq-B-5. dashed: Aq-A-1 
(DM-only). 
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Figure 18. Pair production rates for Aquarius runs re-simulated with baryons (see [53]). solid 
(from top to bottom): Aq-A-5, Aq-E-5, Aq-C-5, Aq-D-5, Aq-F-5, Aq-B-5. dashed: Aq-A-1 (DM- 
only). All lines assume m<f, = 1 GeV. 



4.1 Comparison with Previous Results 

As a final note, we should make a direct comparison with the similar analysis of [22]. 
In their study of the same problem, they concluded negatively that upscattering could 
provide the necessary rate to explain the INTEGRAL 511 keV signal absent the presence 



- 18 - 



A = 1 MeV 



A = 2 MeV 




»V (GeV) ntf (GeV) 

Figure 19. Pair production rates for Aquarius runs re-simulated with baryons (see [53]). solid 
(from top to bottom): Aq-A-5, Aq-E-5, Aq-C-5, Aq-D-5, Aq-F-5, Aq-B-5. dashed: Aq-A-1 (DM- 
only). All lines assume m x = 1 TeV. 

of populated metastable states. We propose two reasons why our analysis reached different 
conclusions from theirs, even when using the same Einasto profile parameters. 

First, they scanned a much broader range of parameters than we did but their results 
are quoted in terms of their dimensionless parameters. While experience has taught us to 
look for resonance regions in scans of and m x , it is not immediately apparent what 
combinations of the dimensionless parameters will yield the same behaviors. Moreover, 
even having identified those, it is not obvious how fine a graining is required to find them. 
For these reasons we believe their scans may have simply missed the specific combinations 
of parameters we study and in particular resonance regions. 

We can compare directly with their results by noting that T = a d m x /A, rj = ad m^j A 
and 28 M = A. Perhaps our best candidate point from the Aquarius A-l profile is = 
1 GeV, m x ~ 600 GeV and A = 1 MeV (recall that throughout this work we have set 
ad = 1/100). Translating this into their variables gives V = 60, rj = 10 and 8M = 0.5 MeV. 
From the top-left group of their figure 5 they only have comparable values for T = 10 and 
r = 100. The two resonance regions we see along A = 1 MeV on the = 1 GeV plot in 
figure 7 would fall between these T = 10 and T = 100 plots. The most direct comparison 
between our result and theirs might be for ms = 1 GeV, m x = 1000 GeV and A = 1 MeV. 
From figure 7 we can calculate log 10 (R/R o b s ) — —0.7. This appears to be compatible with 
their equivalent point ofr = 100, t/ = 10 and 8M = 0.5 MeV. 

Second, we believe the expression for the escape velocity used by [22] is likely a signif- 
icant underestimation. The expression used by [22] is taken from the earlier [19], and has 
a local escape velocity of roughly 400km/s, lower than most current estimates [59]. From 
figure 3.3.3 one can see that our rates drop by a factor of 5-10 for low mass WIMPs and 
2-5 for higher mass WIMPs when we assume this local escape velocity. We believe that it 
is a combination of these two factors that lead us to different conclusions. 

More broadly, our inclusion of the recent results involving baryonic simulations further 
demonstrates that astrophysical effects can naturally boost the rates into the expected 
levels. This would likely affect the qualitative conclusions of [22] as well, but is outside the 
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scope of their work and relies upon results that appeared subsequent to their paper. 
5 Conclusions 

Radiation from electron-positron annihilation in the center of the galaxy has been seen 
since the early 1970's: first by balloon-bourne experiments [8-10] and later by satellites 
such as INTEGRAL. Today we understand there is a clear bulge plus disk morphology with 
the bulge having a full- width half-max of about 8°. The flux from the bulge component is 
estimated to be roughly 10 43 pairs/s though it is unclear how much the disk is contributing 
to the bulge flux. There are many proposed astrophysical sources in the center of the 
galaxy — pulsars, supernovae, LMXBs and microquasars to name a few — but no single 
source seems to have the right flux and shape. All of these sources (except possibly LMXBs) 
are expected to have disk-like morphologies and total flux contributions on the order of 10 43 
pairs/s. XDM scenarios naturally have a bulge-like shape due to the radial dependence of 
the number density. 

In this work, we have performed a thorough numerical investigation into the plausi- 
bility of XDM scenarios explaining this bulge-shaped signal. To do this we first found 
the upscattering cross sections by solving the Schrodinger equation in the basis of partial 
waves for two two-particle states coupled through a Yukawa-type force. After integrat- 
ing the cross sections over the relative velocity distribution to get the thermalized cross 
sections, we then find the rates of pair production by assuming an Einasto profile for the 
WIMP number density. Due to numerical instabilities we were only able to calculate the 
first seven partial wave modes, but we find that in the range of reasonable model parame- 
ters this is sufficient. We expect uncertainties in the number density profile and the local 
galactic escape velocity to eclipse the uncertainty from higher Z-mode contributions. 

Uncertainties in the calculation of electron-positron production rates in XDM scenar- 
ios come from both model parameters and astrophysical parameters. While varying the 
Einasto profile parameters individually leads to monotonic changes in the rates, variations 
of the WIMP mass and force carrier mass are not so simple to predict. We find that using 
the Aquarius A-l Einasto values our pair production rates are on the order of 10 41 -10 42 
pairs/s. These rates could easily change by an order of magnitude or more via minor 
changes to the WIMP mass and profile parameters. When using the Einasto parameters 
from Aquarius simulations including baryons we find rates that are generically of the order 
10 42 -10 pairs/s and even as high as 10 pairs/s. In light of this and the uncertainty in 
the bulge flux contribution from the disk, we believe XDM scenarios can provide a natural 
explanation of the anomalous electron-positron annihilation signal. 
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